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ABSTRACT 


Methods are presented that can be used to make multiple, overset grids 
communicate in a conservative manner. The methods are developed for use with 
Chimera overset method using the PEGSUS code and the OVERFLOW solver. 
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CHAPTER 1 


Introduction 

The field of Computation Fluid Dynamics (CFD), is the study of fluid flow using 
numerical methods to solve equations that govern the physics of these fluids. Tradi- 
tional methods for understanding fluid dynamics have been the use of experimental and 
theoretical methods. However, since the invention of the digital computer, and more 
recently the high-speed digital computer, the field of CFD has grown tremendously in 
both its capabilities and numerical methods. Other factors that have contributed to the 
growth of CFD methods are relative computer costs, increased performance and avail 
ability of computers, including access to high-speed multiprocessor supercomputers all 
the way down to the workstations that have become a common desktop tool to the 
researcher [2,3]. 

Today, two major grid systems exist for solving complex three-dimensional flow- 
fields [4], The difference lies in the way the flowfield is discretized. In the first method, 
the flowfield is discretized into triangular-shaped elements in two dimensions, and tet- 
rahedral elements for three-dimensional fields (Figure 1.1). This type of grid is called an 
“unstructured grid” since the grid points cannot be associated with grid lines. For the 
second method, the flowfield is discretized into quadrilateral elements in two dimen- 
sions, and hexahedral elements in three dimensions (Figure 1.2). In this method, the 
grids are called “structured grids” since the grid points can be associated with grid lines 
in an ordered manner [4] . 

Structured grids can be divided further into two more groups. For complex geome- 
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tries, it is unrealistic to create a single grid. These geometries are usually modeled using 
several grids that adequately resolve specific areas on the geometry of interest. The 
resulting grid system is then combined using some method for communication between 
zones. This is where the two submethods differ. One method is the chimera grid 
approach, or overset method, where grids overlap each other. The second method is 
where grids abut against each other, and are called patched grids or blocked grids. 

For the research presented here, the manner in which patched, structured grids com- 
municate and pass information between each other will be investigated. Furthermore, 
this method will be incorporated into a flow-solver that was developed for overlapping 
grids. There are several reasons for doing so. The method in which overlapped grids 
communicate with each other is strictly by interpolation. The physics of the flow across 
the interface are not computed and thus conservation is not maintained [14]. Under cer- 
tain circumstances, the method can be conservative. However, this is rarely achieved 
when modeling real world complex geometries. Secondly, incorporating a patched grid 
interface will allow users to model regions of the flowfield with the option of a conserva- 
tive interface. This will allow flow discontinuities to pass across the patched interface 
smoothly as if the interface is transparent to the flowsolver. Lastly, it is sometimes con- 
venient to model portions of the flowfield with a patched grid instead of overlapping 
grids that require a certain amount of overlapping to work correctly [14]. This also 
requires that portions of the overlapped grids be cut out of the geometry being modeled 
(Figure 1.3), which is a process that can be extremely time consuming [15]. 

This study will concentrate its efforts in providing a patched grid interface for the 
OVERFLOW flow-solver with the option of making the developed scheme fully conser- 
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vative. OVERFLOW is an in-house code of NASA Ames Research Center for solving the 
thin layer Navier-Stokes equations with the option of using overset grids [12] . OVER- 
FLOW requires the use of a preprocessor program called Pegsus. Pegsus provides 
OVERFLOW the appropriate information on how a set of overlapping grids, or zones, 
modeling a geometry communicate with each other [11]. Pegsus was written for overset 
grids and will need to be modified and new procedures implemented in order for it to 
handle patched grids. 

The current research project was divided into two areas. The first portion of research 
was to investigate how to convert patched grids into overset grids. The method 
employed here is called volume extension and was briefly examined and determined to 
be too time consuming for the end user. Not only does the user have to extend the 
patched grids to create an adequate overlapping region, but also has to use Pegsus to 
create interpolation stencils and hole cuts if necessary. Although the results were some- 
what promising, the volume extension program developed was determined not to work 
under certain circumstances due to grid topology and that an excessive amount of user 
input to the program was necessary. 

The second portion, and the bulk of this research, was to modify and implement new 
methods into both Pegsus and OVERFLOW in order for these programs to handle 
patched grids. This was further divided into two separate development steps. The first 
step was to enhance Pegsus and ensure that the interpolation stencils created for patched 
grids were valid and would work correctly in OVERFLOW. To quickly test this new ver- 
sion without making any significant changes to OVERFLOW, the flow-through bound- 
ary condition was coded. This boundary condition was easily incorporated into 
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OVERFLOW since the method used to update boundary points for overset grids is simi- 
lar to how the boundary points for patched grids with the flow-through condition will 
be updated. The difference lies in where the flowfield information is obtained. For over- 
set grids, a boundary point, or a recipient point, lies within a grid cell of another grid, or 
a donor grid. The node values of this donor cell are used to update flowfield quantities 
strictly by interpolation. For the flow-through condition, the boundary point lies on a 
common plane or interface with the donor grid. Half the update for this recipient point 
will be from the grid it belongs to. The remaining half of the update will be from the 
donor grid using the interpolation stencils from Pegsus. The only requirement is that the 
update information can only come from interior node points of both the recipient and 
donor grids. These enhancements to Pegsus for allowing patched grids and the flow- 
through condition will also be useful when incorporating the conservative patched grid 
interfacing into OVERFLOW. 

The second step for the patched grid capability will be to implement a finite differ- 
ence scheme for patched interfaces into OVERFLOW. Here, one-half of the flux differ- 
encing is computed for each recipient grid. The remaining half, from the donor grid, is 
computed and combined using communication information from Pegsus. Although the 
scheme to be presented is fully conservative for patched boundaries with point-to-point 
matching, it does have further potential of being fully conservative for boundaries with- 
out point-to-point matching. 

The test case for all three methods investigated, the volume extension, flow-through 
boundary condition, and the fully conservative finite difference approach, will be the 
ONERA M6 wing. Comparisons will be made to a single zone case as well as to experi- 
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mental results. 


Page 5 



CHAPTER 2 


Governing Equations of Fluid Mechanics 

The fundamental equations of fluid mechanics are based on three laws, the Conser- 
vation of Mass, Conservation of Momentum, and the Conservation of Energy laws. 
Applying these laws to a fluid flow results in five, coupled, non-linear, partial differen- 
tial equations. The first equation is derived from the conservation of mass and is known 
as the continuity equation. The Momentum equations are derived from the conservation 
of momentum law and represent the momentum components in three orthogonal direc- 
tions. The energy equation is derived from the conservation of energy, which is the First 
Law of Thermodynamics [1]. 

For a cartesian coordinate system, the continuity equation can be written in conserva- 
tion-law form as (Appendix A) [1]: 

dp + 9(p u) + d(pv) + d(pw) = Q (2 j) 

dt dx dy dz 

Here, u, v, and w represent the three velocity components in the x, y, and z directions 
respectively. Lastly, p is the density. The next set of equations are the momentum 
equations, also known as the Navier-Stokes equations. However, it is a common practice 
to include the continuity equation and the energy equation (defined later) as the set of 
equations known as the Navier-Stokes equations [2] . 
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Written in conservation law form, the momentum equations are (Appendix B) [1]: 


d(p u) 3(p u + p) , d(puv) , 3(pww)_ ° X xx , ° X xy , ° X xz , 0 o\ 

“5r + — 51 — + _ 5T + — & 17 + 3y + 3z 1 ' 

3(pv) . d(pH v) . 3(pv 2 + p) 5(pvw) '-^xy 9t )7 p o, 

“5 r + _ ST + 5y + “5^“ 17 + 17 <k 

3(pw) , d(puiv) 3(pvw) 3(pw +p)_ ^xz ^ yz ^ x zz ,,,, 
“5T + “57“ + “5^ + Tz 1T + dy + dz 

Where p is the pressure and x xx , X yy , X zz , x xy , x xz , x yz are the viscous stress tensors. They 
are defined by the following equations: 


_ 2 fr. du dv dw 
X xx 3^\ dx dy dz 


2 f^dv du dw 
3^\ dy dx dz 


2 frydw du dv 
3^V dz 3* dy 


(du dv 


(dw du 

nj; + 5i 


(dv dw 

V - Ka + 3y 
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The last equation, the energy equation, is defined as follows (Appendix C) [1]: 


dE. d((E. + p)u) d((£, + p)v) d((E+p)w) 

1 -| - -|- + = 

dt dx 3y dz 

d ( uX xx + vx xy + WX xz - , d( - UX xy + vx yy + wx yz ‘V, 

dy 

d ( UT xz + vx yz + WT zz-^ ) 

dz 




(2.15) 


where c p is the specific heat at constant pressure, p is the coefficient of viscosity. Lastly, 

Pr is the Prandtl numbers. The Prandtl number is a ratio of the energy dissipated by 
friction to the energy transported by thermal conduction. This can be seen as the ability 
of a fluid to diffuse momentum and energy on a molecular level. 
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Vectorization of Navier-Stokes Equations 


Equations (2.1), (2.2), (2.3), (2.4), and (2.11) can easily be put into vector form. Rewrit- 
ten, they are: 


0£ 9(pn) 9(py) d(pw) Q 

at a* d y d z 

9(p«) 0(pM~ + p ) 3(pMv) + 0(pMVC) _ dT X x + *xy + z 

dt dx dy dz dx dy dz 

d( pv) dipuv) 9(p V 2 + p) 9(pvn0 = <^xy <^yy <^yz 
dt dx dy dz dx dy dz 


3(pw) 0(pMw) 9(pvw) 9(p W^+_p) _ *yz ^ZZ 

dt + dx dy dz dx dy dz 


(2.16) 


dE t d((E t + p)u) 3((E, + p)v) 3((E, + p)vy) _ 
dt + dx + 3y dz 


d(uX xx + vX xy + wX xz -q x ) _ d(ux xy + vx yy + wx yz -q y ) _ 

S3? ay 


^(“ T « + VT y g + WT «-g z ) 

It follows then that the Navier-Stokes equation in conservation-law form may now be 
written in vector form as: 


dE dF„ dG, 


dQ dE dF dG_ v v 

dt + dx + 3y + dx dy dz 


(2.17) 
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Where Q, E, F, G represent the inviscid flux vectors: 


And E 


p 


pu 


pv 


pw 

pu 


2 

pu + p 


puv 


puw 

pv 

E = 

puv 

F = 

pv 2 + p 

G = 

pvw 

2 

pw 


puw 


p vw 


pw + p 

A 


(E { + p)u 


i {E t + p)v 


(E ( + p)w 


F , G v represent the viscous flux vectors: 


0 


0 



V 

V 

< 

II 

V 

^ xz 



uX xx + vX xy + w% xz - *x_ 


“ % xy + n yy + wx yz-1y. 


0 


G, = 


xz 


yz 


zz 


uX xz + vX yz + wX zz-^z 


The next step will be to non-dimensionalize these equations. 


(2.18) 


(2.19) 
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Nondimensional Form Of Navier-stokes Equations 


The Navier-Stokes equations are often put in nondimensional form. There are 
several advantages to doing so. First, characteristic parameters such as the Reynolds 
number, Prandtl number, and Mach number can be varied independently from each 
other. The second reason is that flow variables such as velocity are normalized (the 
resulting values are often between 0 and 1) [3] . For this analysis we will nondimension- 
alize the governing equations with freestream values for the flow variables and a refer- 
ence length for the spatial variables. The following non-dimensionalizing variables will 
be used [3] : 



Where the terms with asterisks are nondimensional. Applying these nondimensional 
variables to Equation (2.17) we get (Appendix D) [3]: 

dQ* dE '* dF ’* + 3G* _ d£* v + ^ F *v + (2 

dt* + dx* + dy* + dz* 3^* dy* dz* 
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Where Q*,E*, F*, G* are: 


Q* 


And E* v , 


p* 

P*M* 

p*v* 

E* = 

p*w* 

P + P* 

p*w* V ; * 

F* = 

P*V* 
p* M* y* 

p * v * 2 + 

G* = 

P* VV* 

p *u*w* 
p* v* w* 

p *w* 


p *W* W* 


P* y* w* 


2 

P*VV*“ + p* 

E* t 


(E* t + /?*)v* 


(E*, + p*)v* 


(E* t + p*)w* 


F* v , G* v . are: 


£* = 


0 
x* 


XX 


xy 


xz 


u* X* + V* X* rv + w*X* X7 - q* 


F* = 


0 


JCJ 

x* 

X yy 


yz 


xz 


«* X* rv + v* X* vv + w* x* yz - q * 


xy 


G* = 


yy 

0 

x* 


xz 


yz 


zz 


U*T* XZ + v*X* V7 + w*X* 77 - q* 


yz 


zz 


( 2 . 22 ) 


( 2 . 23 ) 
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Also: 


E* = p*le* 


u* + V* + w* 


The viscous stress tensors are also given by: 


* _ 2 )J* dw* 

xx IRe^y dx* dy* dz* 


_ 2 n* /^ 3v* du* 

3 Re^V dy * dz *) 


yy 


r * _ 2 jii* /^ 3w* du* _ 3v* N \ 
zz 3 Re dz * 3jc* dy*) 


xy 


- I 11 * f^ u * + d v * ^| 

R e l\dy* dx*) 


= v*_(dw^ 

* z Re L \dx* dz*) 


* = ("gv! , 

y z Rej\dz* dy* ) 


The heat flux terms are given by: 


«*x 




dT* 


(y- \ )M 2 00 Re L Pr dx * 


q* v = - 


dr* 


(y- \)M 1 00 Re l Pr 


dy * 


(2.24) 

(2.25) 

(2.26) 

(2.27) 

(2.28) 

(2.29) 

(2.30) 

(2.31) 

(2.32) 
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We can now rewrite Equation (2.22) in its final form. For convenience we will drop the 
asterisks for the remainder of this work. The Reynolds term has been factored out from 
previous derivations. 


dQ dE dF dG = _L( d _^Z + d Zl + 
dt + dx + dy + dz Re L { dx dy dz ) 


The next step in preparation for applying a numerical scheme to the Navier-Stokes 
equations will be to apply a coordinate transformation to the governing equations. 



Coordinate Transformation 


The governing equations derived thus far have been expressed in a cartesian coordi- 
nate system (x, y, z). This system is called the physical domain. For finite difference 
methods, this type of coordinate system is most efficient. However, for real world geom- 
etries, such as a complete aircraft, an orthogonal coordinate system is unrealistic. These 
type of configurations require a nonorthogonal coordinate system. Therefore it is neces- 
sary to transform the physical domain, (x, y, z) coordinates, to a computational domain, 
(£ , T| , £) [2]. This transformation not only creates and allows an equally-spaced compu- 
tational domain, but also allows the user to align one of the computational directions 
along a specific physical boundary such as an aircraft surface. This will help in applying 
boundary conditions. For this analysis, the Xi direction^ ) is the streamwise direction, 
Eta(t|) is the spanwise direction, and finally, Zeta(£) is the normal or viscous direction. 
Figure 2.1 illustrates the physical and computational domains. For a completely general- 
ized transformation, consider the following transformation: 


% = % z) 

(2.35) 

r| = Ti(x,y,z) 

(2.36) 

N? 

>< 

II 

(2.37) 

X = t 

(2.38) 
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The inverse of these equations are: 


X = x(£>, rj, o 

(2.39) 

y = y(S> 0 

(2.40) 

Z = 0 

(2.41) 

t = X 

(2.42) 

Applying the chain rule of partial differentiation, the partial derivatives of Equations 

(2.39) through (2.42) are: 


a _ a§ a ari_a_ a$_a^ 

ax dxdt, + axati + axac, 

(2.43) 

a _ a^ a ari a ag a 

dy ~ dydt, + aydri + 

(2.44) 

a a^a a^a a^a 

dz dzd^ + dzdn\ dzdC, 

(2.45) 

a a^a a^a a; a «h a 

dt dt a^ + dt dr) + dtd^ dtdx 

(2.46) 

The metrics in these equations are: 


^ „ ^5 _ r 

dx ~ s *’3x '*’3x s * 

(2.47) 

^ _ % ^1 _ ry 55-r 
ay ” ^y’dy 'y’dy s y 

(2.48) 

dz " dz ' Z ’dz Sz 

(2.49) 

dt ~ ’ r 3r " ^<’3t ^’d» ' 

(2.50) 
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These metrics, in general, cannot be determined using analytical methods. Therefore, a 


numerical method must be used. First, consider the following differential expressions: 


dt = t T dx + t^d^ + t r[ dy\ + t^dC, (2.51) 

dx = x x dx + x^dL, + x^dx\ + x^dC, (2.52) 

dy = y x dx + + y^dj] + y^dC, (2.53) 

dz = z % dx + z^dL, + z^dr\ + z^dC, (2.54) 


From Equation (2.42), time, t, is only a function of the computational domains’ time, t . 
Therefore the following partial derivatives must equal zero. 


3 ? _ dt_ _ 3 /. 
a^ " ari “ ac 


(2.55) 


Also using the relationship in Equation (2.42), the partial derivative of t with respect to t 
must equal one. With this information and Equation (2.55), we can rewrite Equation 
(2.51) as: 

dt = dx (2.56) 

Placing Equations (2.51), (2.52), (2.53), and (2.56) in matrix form we have: 



(2.57) 
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Reversing the role of the independent variables, it can be shown the reverse is true with 
the following differential expression: 


dx = dt 

(2.58) 

d\ - % ( dt + \ x dx + ^ y dy + \dz 

(2.59) 

dy = r\ f dx + + r\ y dx] + x\d ^ 

(2.60) 

dz = + S/H + 

(2.61) 


In matrix form we have: 



(2.62) 


From Equations (2.57) and (2.62), the following must be true: 

10 0 0 

Sy 

1/ ^ ’ly \ 

t, ^ 

Using a symbolic math program such as Maple V, the right-hand-side of Equation (2.63) 
was solved for (Appendix E) [3]. The transformation metrics can now be solved for with 
the following results: 



10 0 0 


't x £, x n x ? 


^ ^ 


^ 


-i 


(2.63) 


= ; (Vr :y { z n ) 

ly = *■<{■$ 
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(2.66) 


'Hjp = J(y^-y^z^) (2.67) 

T|^ = J(x^z^-x^) (2.68) 

r\ z = J(x^y^-x^y^) (2.69) 

t, = (270) 

?, = • / (vs _ - t 5V (2 - 71) 

C z = J(x^y^ - x^y^) (2.72) 

= -(^ + ^, + 7^) < 2J3 > 

ii=-(¥i + )'A +I tV (2J4) 

= + < 2 - 75 > 


Equations (2.73), (2.74), and (2.75) can be rewritten using the previous spacial metrics as: 

\ t = ^(j,v 7 i;V t ^'Vc )tl i(Vr , i; ,, ii )1 (276) 

t| ( = -y[.* x (y^-E - y^Zf-) + y x (x^ - zcj-z^) + z t (jt^ - xyyrj 1 (2 - 77) 

C, = - y^) + y t (j: n z 5 - *j=z n ) + z t (^y„ - V$> ] (2J8) 

Where the Jacobian, J, is defined as: 

j _ giizILO = I (2.79) 

d(x ,y,z) x ^(yn z ^~y^ z ri^ + ^~ x ri^ z ^~ y ^ z ^ + x ^ y ^ z r\~ y ri z ^ 
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Applying this generalized transformation to the governing equations we get (Appendix 
D) [3]: 


dQ dE 3F 3G_ 1 ( dE y aG v 

3x + 3^ + 3r| + 3 C, Re L 3£, 3rj 3^ 


(2.80) 


"pi p" 

p u P uU + i > X P 

ip v E = j pvU + l y p 
pw p wU + ^ z p 

Al (E l + p)U-\ t p 


pV 

puV+T] x P 


pw 

puW + t, x p 


pvV + rj^p G = J P vW + ^ y P 

pwV + r\p pwW + C,p 


{E t + p)V -r\ t p 


pwW + C, z p 
( E t + p)W-C >t p 


(2.81) 


Where U,V, and W are the contravariant velocity components defined as: 

u = §, + 1,? + 5/ + ^ v 

V = 11, + V + V + V (282 ) 

w = C, + + ; z v 

The contravariant velocities represent the velocity components that are perpendicular to 


planes of constant ^ , r\ , £ [3] . 
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And the viscous flux terms E v , F v , G v are: 


K = 

K = jClA + l/v + le 0 .) 

g v = j«A + V» + W 

Where the viscous shear stress terms are given as: 

■ i3fe- [2(? A +1, A + t A ) - (? » v 5 + Vti + Vc ) - 

L 

^z w % + \ w n + ^z w 0 ] 

tyy = l^ [2( Vs + Vn + Vc M Vs + V n + W- 

L 

(^ + Vi + W 1 

= 5^ 12( V + VV V>-<V + VV 

L 

( Vs + Vn + i; A )1 

V = ;?r ( Vs + Vn + V? + Vs + 11 * v n + V;* 

Zi 

T « = j£r< Vs + W Vs + Vs + 11 *% + Vc> 

L 

v = «r ( Vs + Vn + Vs + Vs + Vn + Vs> 

L 


(2.83) 

(2.84) 

(2.85) 

( 2 . 86 ) 

(2.87) 

(2.88) 
(2.89) 
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Also, the heat flux terms are given by: 


= ~ 


- 2 

(Y-l )M „Re L Pr 


q y = 


1 + W ?/;) 

(y- i)M ooRejPr 


1 


— \ (^h + \ r n + V0 

(y- 1)M ooRe L Pr 


(2.90) 


(2.91) 

(2.92) 
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Thin Laver Approximation 

In order to adequately resolve viscous gradients due to solid surfaces, many grid 
points must be clustered in these regions. For flows in which separation is minimal, gra- 
dients of the stress terms normal to viscous surfaces have been found to be much larger 
than gradients parallel to these surfaces [10]. Therefore, it is unnecessary to construct 
grids with fine resolutions in directions parallel to viscous surfaces. To expand on this 
analysis we must include computer requirements and limitations. Since a large amount 
of computer time and storage is used in resolving these gradients, any valid reduction in 
the full Navier- Stokes equations would greatly increase the efficiency and productivity 
of the researcher. It has been shown that if the grid spacing parallel to viscous surfaces is 
too coarse, and that the full Navier-Stokes equations are being solved, the resulting solu- 
tion indicates that viscous gradients in these directions have not been fully resolved 
[2,10]. It then only makes sense to drop the viscous terms in which partial derivatives 
have been take with respect to directions parallel to viscous surfaces. As mentioned 
before, after the governing equations have been transformed, it is a common practice to 
have the viscous direction to be the Zeta direction. Therefore, any of the viscous flux 
terms that contain partial derivatives with respect to Xi or Eta can be dropped since these 
terms have been shown to be much smaller than those terms that have partial deriva- 
tives with respect to Zeta. Returning to Equation (2.80), we can now drastically reduce 
the viscous flux terms by dropping all terms with derivatives take in respect to Xi and 
Eta. The resulting equation is: 



Where s is defined as (Appendix F): 


5 = 7 


+ Cv + (£)“C + + C y v c + 5^5) 

n(C, 2 + Cy + 5*)v ? + + c y v 5 + c^ t ) 

+ 5^z(C J ,» ? + 5 y v; + ^z w t) 

r 2 r 2 r 2T fdw 2 dv 2 dvv 2> \ dr"! 1 f r 9 m r dv r dw\ r , r r v 

(C ( + C, + ^Kac + ac + irJ + “ac J 5^55 + C '^ + ^r* u + C ' v + ^ w) 


(2.94) 
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CHAPTER 3 


Volume Extension 

The idea behind volume extension is to extend patched grids so as to create an over- 
lapping region between grid zones. The overlapping region is required in order for Peg- 
sus as well as OVERFLOW to work correctly. There are several reason for investigating 
this approach. First, if this method of creating an overlapping region works and with lit- 
tle user input, it would alleviate any need to modify PEGSUS and OVERFLOW. This 
would be the major reason. Secondly, many researchers create patched grids for use in 
their own flow-solvers. Often, these researchers require validation of their results using 
other flow-solvers, such as OVERFLOW, that use the overset grid technology. Instead of 
completely remodelling the configuration, the existing patched grids could be trans- 
formed into overset grids. As will be shown, the program developed, “patched-2-over- 
set”, was not as successful as hoped. 

At this point it is a good idea to explain how grids are setup and how the relationship 
between the physical plane and computational plane are determined. As previously dis- 
cussed, the physical plane is in x, y, and z coordinates in space and the computational 
space is in Xi, Eta, and Zeta coordinates. These Xi, Eta, and Zeta coordinates are most 
often indexed using J, K, and L indices. Constant Xi planes all have the same J index. 
The same is true for the Eta and Zeta planes both having constant K and L indices. The 
indices are integers starting at one and running to however many grid points are in that 
specific direction. Figure 3.1 demonstrates this idea. 

For the program to work successfully, the user needs to identify several features of 
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the grid to be extended. But first, some assumptions will be made. The program 
assumes that the boundaries of the interface will match up perfectly with the donor side 
patched grid boundary. The computational indices do not necessarily have to be the 
same, but the physical boundaries of the interface need to match up. This assumption 
was made since the test case of the ONERA M6 wing had this feature. In fact, this test 
case was actually a single grid as shown in Figure 3.2. It was then chopped into four sep- 
arate zones, thus creating twelve, mesh continuous, patched interfaces as shown in Fig- 
ure 3.3. However, this assumption would rarely exist for complex geometries where 
patched interfaces may have several grids as donor interfaces. 

The user then needs to provide information to the program in order for it to under- 
stand where in the computation space the interface exist. This includes which grid the 
interface exist in. Also, the J, K, and L index ranges that define where in the computa- 
tional space the interface lies. As one will assume, the interface will always be either on 
a J=l, J=Jmax, K=l, K=Kmax, L=l, or an L=Lmax plane. The next input is tell the pro- 
gram where the donor side interface lies. This will requires providing the same infor- 
mation as previously given to the program for the recipient side interface. To help 
reduce user input, it was assumed that the reciprocal is true. When the donor side inter- 
face is provided, it is assumed that it will be a future recipient interface with its previous 
recipient interface now becoming the donor interface. 

Providing both interfaces is very important, the program uses the donor grid inter- 
face boundaries as a stencil on how to extend the recipient grid. It also uses the plane of 
points behind the donor interface as a reference plane on how far to extend the recipient 
grid and how the recipient grid will conform to the donor grids outer boundaries. Fig- 
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ure 3.3 demonstrates this idea in two-dimensions. However, it will be shown that this 
may cause problems under certain circumstances. The user also has the option of how 
many grid points to extend the grid. The grid is extended using the donor grids inter- 
face plane and all subsequent planes behind the donor interface depending on how 
many planes are added. Since the recipient grid is extended using the donor grid, the 
user specifies a maximum stretching ratio, in returns the program alerts the user if the 
stretching ratio is exceeded. Since there should be little stretching across the interface, 
this alert should help the user in problem areas of the grid to be later refined. Table 3.1 
summarizes the user required inputs as the program is executed. 

The test case chosen for this research is the ONERA M6 wing as previously shown in 
Figure 3.3. Table 3.2 details the grid system information. For this test case, we will dem- 
onstrate the use of the extension program on one of the patched surfaces. It will be 
shown that all twelve interfaces will have similar results. For this demonstration, the K 
maximum patched surface on zone three will be extended one grid plane using the 
donor grid for interpolation. As Figure 3.4 shows, the K maximum interface is also on 
the same physical plane as another patched interface of zone three, the J minimum sur- 
face. The same is true for the donor grid, zone two. The K maximum surface of the 
recipient grid matches up with the K maximum surface of zone two. Also, the J mini- 
mum surface of zone three matches up with the J maximum surface of zone two. Let’s 
now concentrate on the K maximum surface of zone three, our demonstration surface. 
However, all four surfaces just mentioned will in one way or another affect this surface, 
and will lead up to one of the setbacks of volume extension. Instead of trying to visual- 
ize this extension in the physical plane, it will be much easier to show in the computa- 
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tional space. Referring to Figure 3.5, the K maximum surface and the J minimum 
surfaces are shown. This surface will now be extended as shown in Figure 3.6. The 
extended plane is shown in red. This will now add one more dimensional grid point for 
the K direction. Disregard for now using the donor grid for the distance and placement 
of this extended surface. Remember that this was a single grid and was split at symme- 
try planes, so using the recipient grid for an extension direction and distance will have 
the same effect as using the donor grid. Now, lets extend the J minimum surface without 
extending the new K maximum plane. Figure 3.7 shows the results to this with the new 
J minimum plane in blue. As a result of extending the K maximum, an extra point has 
been added. As shown in Figure 3.7, the green dots are the result of extending the first K 
maximum plane and the second J minimum plane. The end result is a series of new cells 
at the corner of the computational grid. This is fine for the computational space since the 
angle alpha between the new K maximum plane and the new J minimum plane is ninety 
degrees. However, referring back to Figure 3.5, both the K maximum plane and J mini- 
mum plane are all on the same physical plane. This will make the angle alpha zero 
degrees, and as a result the new set of grid cells at the corner of the computational plane 
will collapse on themselves in the physical plane. An undesirable result that creates zero 
and negative cell volumes. 

It just happens that all twelve patched surfaces in the test case have this feature. 
There are several options to further investigate. The other option available in this pro- 
gram is to use the recipient grid to extend the patched interface. That is use the plane 
behind the patched interface and the interface itself to linearly extrapolate a plane. As 
stated earlier, the four zones were created at symmetry planes and as a result, the two 
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options created in the program will have the same effect. Another option is to extend the 
interface normal to the patched surface. This again will have the same effect as before 
with corner point B being coincident with points A and C as shown in Figure 3.7. The 
last option used was to go ahead with the original grid extension and force the angle 
alpha to be something other than zero degrees. This will force all the points beneath 
points A and C to spread outward from the points beneath point B. In the end this 
worked fine and the grids were run through Pegsus with no problems. 

This test case is rare in that all the patched interfaces had this feature. For the most 
part, this program could handle most real world complex geometries modeled using 
patched grids. However, most often, a patched interface will not be an entire computa- 
tional plane. Usually the patched surface will be a subset of a computational plane with 
boundary conditions applied to the remaining points on the surface. This will create 
problems later on for the user once the entire plane is extended. The extended portion of 
the plane that had boundary conditions will now have to be cut out during the Pegsus 
procedure. This may be a terribly time consuming task if for say the test case has eighty 
patched grids with a good portion of them needing hole cuts due to the grid extension. 
Also to consider, even if extending patched grids to create an overset grid system will 
still leave the interfaces unconservative. Due to the outcome of this investigation, and 
the reality of having further problems with grid extension, it was decided to continue in 
another direction. The results for this test case will be presented later in Chapter Six. 
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CHAPTER 4 


Pegsus Modifications 

As mentioned earlier, in order to model real world complex geometries, multiple 
grids must be created that model specific regions of the flow field accurately. This proce- 
dure of domain decomposition is called the chimera grid approach [11]. Since the main 
purpose of the flow-solver, OVERFLOW in this case, is to numerically solve the govern- 
ing equations, a separate preprocessing program is required to provide information to 
OVERFLOW on how each grid will pass flow-field information between themselves. 

For OVERFLOW, the program Pegsus is used to provide grid communication informa- 
tion for overlapping grids. Pegsus is also used to cut holes in grids in areas that may 
cause interference with other grids such as overlapping solid surfaces. 

The basic procedure for a researcher to apply the chimera grid approach is to first 
model the flow-field using multiple overlapping grids (Figure 1.3). Next, the user pro- 
vides Pegsus with these separate grids and an input file that provides all necessary infor 
mation on which outer boundaries on each grid will require Pegsus to compute valid 
stencils. A valid stencil is a when Pegsus has determined that a point lies within a donor 
grid cell and that the three required interpolation coefficients are between zero and one. 
For grids requiring hole cuts, the user specifies surfaces in one grid that will make a hole 
in another grid. Pegsus will then determine how to cut the hole and which fringe points 
on the hole will require interpolation stencils. Pegsus in returns creates a single multiple 
grid file containing all grids and an interpolation file [11]. 

For this research, the process in which Pegsus goes about searching for valid stencils 
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and creating hole cuts will not be discussed in detail but rather how to modify Pegsus to 
handle patched grids. In fact, hole cutting is not required for patched grids since there 
are not overlapping regions. 

First there is the issue of the interface itself. For grids of varying resolution or non 
point-matched grids, the interface is not exactly a clean smooth curve where there are no 
gaps or overlapping. For most cases there exist smalls gaps in some regions and small 
overlaps in others as demonstrated in Figure 4.1. The interface points that overlap the 
donor grid and those points very close to the donor grid, Pegsus will have little trouble 
finding valid stencils. It is difficult to determine what is a very close point in terms of a 
specific distance. It will vary from run to run and the grid resolutions involved in each 
run. One of the inputs to Pegsus is for the user to determine and supply a parameter, 
epsilon, that will allow Pegsus to label a stencil as valid depending on how far that point 
may lie outside of a donor grid cell [11]. This will result in interpolation coefficients out- 
side the range of zero and one or that the point lies outside the donor cell by some small 
amount. For the points outside the domain of the donor grid, Pegsus will have trouble 
finding stencils. For these points to receive the appropriate interpolation coefficients, 
they must be projected onto the surface of the donor grid. The projection will not be per- 
manent, it will only provide a physical coordinate to compute the interpolation coeffi- 
cients. This projection requirement will be one of the modifications to Pegsus. A search 
routine will also be needed to find out where on the donor grid interface the projection 
will occur. 

With both the projection and search routines added, some other features will be 
added to speed up the search and add functionality for the patched grid option to Peg- 


Page 31 



sus. Table 4.1 summarizes these additions. 

The first entry in Table 4. 1 is the patched grid option. This option, specified in the 
input file provided by the user, enables the new search and projection subroutines. This 
option is only enabled if a point is labelled by Pegsus as an ‘Orphan Point’. An orphan 
point is a point that Pegsus could not find a valid stencil for. Once a point is labelled an 
orphan point, the new search routine begins. 

The reason for the new search routine is to find the donor grid cell that the recipient 
point lies within or has the possibility of being projected onto in order for interpolation 
coefficients to be calculated. The first step of the search is to find the closest point in the 
donor grid to the recipient point. This donor test point will be the starting point of the 
next search. The next search procedure is a clipping test. This test determines if the 
recipient point has the possibility of falling within a two by two cell face area of the 
donor grid with the donor test point just found as the center point. This is achieved by 
first determining the minimum and maximum x, y, and z values of all nine points that 
make up the two by two cell faces of the test region of the donor grid. With these mini- 
mum and maximum x, y, and z values, the clipping test is passed if any of the x, y, or z 
components of the recipient point lie within the minimum and maximum ranges of the 
donor test region. The clipping test is rather fast and reduces the number of regions in 
the donor grid for the next search procedure to test in. The clipping test is demonstrated 
in Figure 4.2. If the test fails, the search continues in a radial direction from the initial 
test point. Once the clipping test is passed, the next test procedure is executed. 

The final search procedure is a triangle test. That is, the recipient point is projected 
onto the plane of the two by two cell face area using the center point as a reference point. 
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Even if the four cell faces are not on one distinguishable plane, the center point will be 
common to all four planes and is used to project the recipient point the proper distance. 
Once the point is projected, a triangle test is implemented. This test consist of calculat- 
ing four triangular areas and making a comparison. Referring to Figure 4.3, point R, the 
recipient point, is projected onto the donor grid interface test region. The next step is to 
calculate the area of triangle abc. This will be the base area. Next using the new coordi- 
nate of point R after the projection, the areas of triangles abR, acR, and cbR are calcu- 
lated. If the sum of these component areas sum up and equal to the base area, the search 
is over and the grid cell in the donor grid has been found. If not, the triangle test contin- 
ues for this test region. If all eight possible triangle tests fail, the search continues on 
with the clipping test again. Once the triangle test has passed, the donor cell and pro- 
jected point R are passed back to the main search of Pegsus and interpolation coefficients 
are calculated. To help speed up future searches of orphan points passed to the projec- 
tion subroutines, the previous stencil is stored. Instead of calculating the closest point in 
the donor grid, the stored stencil is used as a starting point. Calculating the closest point 
in the donor grid is CPU intensive, avoiding having to use it drastically speeds up the 
search. Since orphan points tend to be confined or clustered to specific regions of the 
grid, storing previous stencils most often reduces calling the closest point search to only 
two or three times depending on grid topology. 

The next option in Table 4.1 is the DGRID option. The user is responsible for provid- 
ing Pegsus a set of possible donor grids. For patched grids, usually only a subset of 
these grids are possible donor grids for one interface and the remaining subsets are pos- 
sible donor grids for the other interfaces of the recipient grid. The criteria Pegsus uses 
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for a recipient point to even start a stencil search in a donor grid is if the points lies 
within the min/max box of that donor grid. However, in some cases and especially with 
patched grids, the recipient point may be in the min/max box of the donor grid but may 
not be its true donor grid. Therefore, the DGRID parameter is an option to further subset 
the user provided donor grid list for each of the interfaces on the recipient grid in ques- 
tion. This has shown to improve speed of execution in most cases. This is because once 
Pegsus has started its search and the point will require projection, Pegsus will continue 
the search over a great deal of the donor grid until it labels it an orphan point. This is 
costly and can be avoided if the user uses the DGRID option. 

The third entry in Table 4.1 is the PROGRD option. PROGRD is a separate program 
to Pegsus and is used to project subsets of a recipient grid onto a reference plane of a 
donor grid [13]. This is primarily used in the viscous regions of grids that have varying 
resolution and there exist a mismatch as to were the solid surface is modeled. Referring 
to Figure 4.4, both the red and blue grids accurately model the cylinder. All surface 
points lie on the solid surface of the cylinder, however since the two grids have varying 
resolutions there exist a mismatch in the boundary layer of these two grids. Referring to 
the blue grid and line A, the point on the solid surface is actually within the boundary 
layer of the red grid. If these were two patched interfaces, the flow-field information 
being passed back and forth between the red and blue grids, in the sensitive region of the 
boundary layer, would be incorrect and will cause discontinuities in the flow. To fix this 
problem, PROGRD will project the points in line A on the blue grid onto the solid sur- 
face characterized by the red grid. The PROGRD option, if specified, will read in the 
resulting projected grid from PROGRD as the set of physical coordinates used to find 
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stencils in PEGSUS. However, the original grid is used in the flow-solver since it accu- 
rately models the geometry. 

The fourth entry in Table 4.1 is the METHOD option. This will tell Pegsus were to 
store the weights for the patched grids interface points. This will allow OVERFLOW to 
distinguish between the two methods used to update interface points. As mentioned 
earlier, the two methods used to update the interface points are the flow-through bound- 
ary condition and the cell-vertex finite volume method. Since the flow-through bound- 
ary condition only uses the first interior points behind the donor grid interface for flow- 
field information, it was decided to make the interpolation coefficients indicate that all 
the weight is on the interior plane. This will require making one of the weights zero or 
one. Depending on what computation plane the interface lies on, one of the weights will 
either be zero or one. For example if the interface is on the J maximum plane, the J 
weight will be zero, thus indicating no weight on the J maximum surface. Since the 
interface is two-dimensional, the weight being set to zero or one is in the third direction 
which has no effect on computing the final weights which are dependent only on the 
remaining two weights. The opposite is made for the cell-vertex method. One of the 
three interpolation coefficients will be set to zero or one indicating the remaining two 
weights are on the interface. For the previous example, the J weight will be set to one 
placing full weight on the interface. When OVERFLOW reads in the weights and the 
stencil, it will be able to differentiate between the two methods depending on where the 
weights are. 

The last entry in the Table 4.1 is for the coincidence check. For the example of the 
ONERA M6 wing, the patched grids are mesh continuous. That is, the grid points line 
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up, or are coincident with the donor grid interface points and vice versa. When the coin- 
cidence check is enabled Pegsus will check for coincidence with the donor side test 
point. This is a fairly easy task of comparing physical coordinates and is the first test in 
the search procedure. If the user knows that the patched grids are mesh continuous, 
enabling this option will speed up execution tremendously. 

The modifications to Pegsus were straight forward. The bulk of the work was creat- 
ing the new search and projection subroutines. The remaining options listed in Table 4.1 
were consequences of debugging the program. Although the speed of the program has 
decreased some what as compared to the original version, the resulting number of 
orphan points is minimal, which is one of the goals for a patched interface. 
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CHAPTER 5 


OVERFLOW Modifications 

OVERFLOW is a three-dimensional unsteady compressible thin-layer Navier-Stokes 
flow-solver developed at NASA Ames Research Center. The development of OVER- 
FLOW has been ongoing for many years and begun as a rewrite from two previous 
NASA codes, the ARC3D and F3D codes. OVERFLOW has several options for right- 
hand-side calculations [12]. These include central differencing, flux split in J (central in 
K, L), Liou AUSM flux split scheme, and finally the Roe upwind scheme. Left-hand-side 
options are the ARC3D 3-factor block tridiagonal scheme, F3D two-factor scheme, 
ARC3D 3-factor diagonal scheme, and also the LU-SGS algorithm. Several turbulence 
model options exist also. These include the Baldwin-Lomax, Spalart-Allmaras models 
and the k - e turbulence model. For this research, the right-hand-side calculations will 
be limited to central differencing and the ARC3D 3-factor diagonal scheme for the left- 
hand-side [13]. 

The first modification to OVERFLOW was to include the flow-through boundary 
condition. This was included since it is a simple modification to the program and was 
used to quickly debug problems in the modifications to Pegsus. Although the flow- 
through boundary condition is one solution to including the patch grid capability in 
OVERFLOW, it will prove to be unconservative. The goal is to include a patched grid 
capability with the option of further enhancing this method to be fully conservative. The 
flow-through boundary will not be capable of being modified to be fully conservative. 
The initial results were very promising, proving that the modifications made to Pegsus 
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for patched grids were working. Also, the results of the ONERA M6 wing showed that 
the flow-through boundary condition worked well with very little flow-field discontinu- 
ities at the interface. This test case did not place patches near any known flow-field dis- 
continuities such as shocks. Later test cases, with patched interfaces placed at a shock 
location on the top of the wing, showed that this method would not allow the shock to 
properly pass through the interface. 

The next modification to OVERFLOW is to include a cell-vertex finite volume scheme 
for the interface points. This method has the further capability of being fully conserva- 
tive. With the interior scheme being a finite difference method, including the cell- vertex 
finite volume method took careful planning and consideration to allow both methods to 
work together. Both right-hand-side central differencing and the left-hand-side ADI 
operator were implemented in several steps. This method proved to be unstable for cer- 
tain test cases involving an initial garbage solution. The initial garbage solution tests to 
verify if the garbage input would pass through the interface cleanly and wash out the 
exit plane. The ONERA M6 wing test did not work at all. Severe pressure oscillations 
occurred immediately ahead of the leading edge of the wing. Further investigation of 
the cell-vertex method is needed to stabilize the scheme. 
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Flow-Through Boundary Condition 

The flow-through boundary condition works very similar as the C-grid boundary 
condition [12]. The C-grid boundary condition is applied to a plane of coincident points 
in the region of a C-grid that folds over onto itself. Figure 5.1 demonstrates a C-grid 
modeling a cylinder. The wake region of the C-grid is were the C-grid boundary condi- 
tion is applied. The update of the boundary points using this boundary condition 
requires using flow-field information from interior points on both sides of the C-mesh 
region. The only requirement is that the points at the boundary must be coincident with 
its corresponding point at the other end of the C-mesh region as demonstrated in Figure 
5.2. In Figure 5.2, the C-grid boundary condition is applied to the J values of one to 
twenty-three. It is automatically assumed the corresponding coincident points, denoted 
by a negative sign, will have the same requirement. Figure 5.3 demonstrates how to 
update the boundary points. For the example, point thirteen is updated using the one- 
half the value of the first interior point behind point thirteen, point Q(13,2), plus one-half 
the value of the first interior point behind its corresponding coincident point, point Q(- 
13,2). The same will be true when updating point Q(-13), it will have the same value as 
its corresponding coincident point, point Q(13). All the points in the C-mesh with the 
boundary condition will be updated after all the interior points have been updated. 

To apply this concept to a patched grid is very similar. For the C-grid boundary con- 
dition, the user specifies only one side of the range of points it is to be applied to. The 
corresponding coincident points are easily obtained since they are in the same grid and 
will also take on the boundary condition. For a patched grid, Pegsus will provide all the 
necessary information on how to identify the donor grid information. This will include 
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providing a pointer and three interpolation coefficients. The pointer identifies the donor 
grid number and the J, K, and L index location inside that grid as to what cell will pro- 
vide the donor side flow-field information. The pointer indices are the minimum index 
of the donor grid cell in each direction (Figure 5.4). Figure 5.4 demonstrates how to the 
recipient point is updated. Pegsus provides the pointer and the three weights which can 
be used to calculate all the weights for each of the node values for the cell. The recipient 
side contribution to the update is one-half of Q { . Where Q { is the first interior point 

behind the interface. The donor side information will be one-half of the sum of several 
weighted Q values. As stated earlier, the weights on the interface points, weights one 
and four, are forced to be zero during the Pegsus process in order for OVERFLOW to 
distinguish the method. 

The modification to OVERFLOW was straight forward. Referring to Figure 5.4, if 
point Q r was an overset point, its update would take on the following formula. 

Q r = +(0 2 Q 2 + (0 3 Q 3 + (0 4 Q 4 } ( 5 . 1 ) 

The update using the overset method does not using any flow-field information from the 
recipient grid. For the recipient patched interface points, the donor side information has 
been satisfied already since weights one and four are set to zero. All that is needed is to 
add the recipient side Q value and divide the total update by two. Therefore, as the chi- 
mera subroutines are iterating over all the outer boundary points, the recipient side Q 
value is added and the total Q value is then divided by two. 

To allow for both overset and patched grids was also a straight forward task to incor- 
porate into OVERFLOW. Pegsus was further modified to insure that patched interface 
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points will always have at least one of the interpolation coefficients set to zero or one 
(only two coefficients are needed, interface is two-dimensional). Then, to distinguish 
between the patched interface points, the overset points will always have values 
between zero and one. If by chance Pegsus found a stencil that had an interpolation 
coefficient value of zero or one, a check is made that will add or subtract a small delta 
value. This will have little effect on the solution but will help OVERFLOW distinguish 
between the two methods. 
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Table 3.1: Patched-2-Overset Required Inputs 


Patched-2-Overset Input Questions 

Iterations 

Enter the number of patched grids in system. 

a=N 

Enter the number of patched surfaces for grid (a). 

b=l..M 

Enter the J, K, and L indices for patch (b) of grid (a) 


Type of Interpolation for patch (b) of grid (a): 
l=straight extrapolation follow grid lines of recipient grid. 
2=straight interpolation, follow grid lines of donor grid. 

1 or 2 

Number of additional planes for patch (b) of grid (a). 

1,2,3 

Enter the number of the donor grid. 

c=N’ 

Enter the J, K, and L indices for donor patch (c) for grid (a) and 
recipient patch (b). 



Table 3.2: Four Zone ONERA M6 Wing Summary 


ONERA M6 Wing: Four Zone, Mesh 




Continuous Grid System 




Zone One: 49x25x33 grid points 

Lower downstream 

grid 

Patched Surface 

J=l,48 

K=2,24 

L=33,33 

Patched Surface 

J=l,48 

K=25,25 

L=2,33 

Patched Surface 

J=l,l 

K=2,25 

L=2,33 

Zone Two: 73x33x49 grid points 

Lower wing grid 

Patched Surface 

J=l,l 

K=2,33 

L=2,48 

Patched Surface 

J=73,73 

K=2,33 

L=2,48 

Patched Surface 

J= 1 ,72 

K=33,33 

L=2,48 

Wing Surface 

J=l,73 

K=2,33 

L=l,l 

Zone Three: 73x33x49 grid points 

Upper wing grid 

Patched Surface 

J=l,l 

K=2,33 

L=2,48 

Patched Surface 

J=73,73 

K=2,33 

L=2,48 


Page 44 












Table 3.2: Four Zone ONERA M6 Wing Summary 


ONERA M6 Wing: Four Zone, Mesh 




Continuous Grid System 




Patched Surface 

J=l,72 

K=33,33 

L=2,48 

Wing Surface 

J=l,73 

K=2,33 

L=l,l 

Zone Four: 49x25x33 grid points 

Upper downstream grid 

Patched Surface 

J= 1 ,48 

K=2,24 

L=33,33 

Patched Surface 

J= 1 ,48 

K=l,l 

L=2,33 

Patched Surface 

J=l,l 

K=l,24 

L=2,33 
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Table 4.1: Pegsus Modifications Summary 


Pegsus Modifications 

Description 

Option 

TYPE = ‘Patched’ 
or 

TYPE = ‘Overset’ 

The ‘Patched’ option will enable the new 
search and projection subroutines 

Required 

DGRID=’donor grid #’ 

This will further subset the donor grids. 
Will only search those grids specified 

Optional 

PROGRD = ‘progrd name’ 

This option will read in grids that result 
from the program PROGRD. Stencils now 
will use the physical coordinates from the 
progrd grid for stencil searches, but retain 
the original grid for the flow-solver. 

Optional 

METHOD = ‘Flow-through’ 
or 

METHOD = ‘Cell-vertex’ 

This option will determine where interpola- 
tion weights will be stored. This will help 
OVERFLOW distinguish between the two 
interface methods, the flow-through bound- 
ary condition and the cell-vertex method.. 

Required 

COINCK = ‘TRUE’ 
or 

COINCK = ‘FALSE’ 

This option will enable or disable the coin- 
cident point check option. Program execu- 
tion speeds up with this option enabled for 
mesh continuous grids. 

Optional 
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Figure 1 . 1 Unstructured surface grid on aircraft geometry 



Figure 1 .2 Structured surface grid on aircraft geometry 
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Figure 2. 1 Generalized transformation between the 
physical and computational domains 
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Figure 3.5 Zone three, patched interface at the K maximum surface 

and the J minimum surface 



Figure 3.6 Zone three, extending the K maximum surface 





and J minimum surface 


Page 53 







Figure 4.3 Triangle test example 
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Figure 5.3 C-grid boudary condition formula 


= 205 + 2 {o) 2^2 + (O 3^3J 


= ( 0 4 = 0 


\ g > 2 “ 3 


Pointer . 


Interface 


Recipient 

Grid 


Figure 5. 4 Patched grid interface with flow-through condition 
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APPENDIX A 


Derivation of the Continuity Equation 



Figure A. 1 Differential fluid element with mass flow rate quantities 
prescribed for each of the elements six faces 

The Conservation of Mass law basically states that mass cannot be created nor 

destroyed. We can demonstrate this by creating a control volume of arbritrary size in 

three dimensions as illustrated above in Figure A.l. Then using the conservation of 

mass law, we can say that the rate of change of mass within the control volume must be 

equal to the net rate of inflow of mass, or in equation form we have: 

dm , A ! 

~ m enter m exiting ^ 

Substituting the mass flow rates for each direction, as illustrated in Figure A.l, we 
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then have: 


( pudzdy) enter -(pudydz + ^j^ldydzdx) 

v OX sexit-lxdirection 

\ (pvdxdz) t -( pvdxdz + ^ V -dxdzdy) + 

L V oj Jexitlydirection 

[ (p wdxdy) emer -[ p wdxdy + dxdydz) 

L \ az J exit-izdirection 

d(p )dxdydz volumemass 

dt 


(A.2) 


Further reducing we have: 


3(pw) 

dx 


dxdydz 


j~9(^vv) dxdydz 


-ixdirection 


zdirection 


, dxdydz] 

L dxy -lydirection 

dxdydz] 

- dt J volumemass 


(A.3) 


Finally, with further reduction we get equation (2.1): 

Bp 3(piQ 3(pv) 3(pw) _ 0 
dt 3x 3y dz 


(A. 4) 
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APPENDIX B 


Derivation of the Momentum Equations 



Figure B. 1 Differential fluid elements with force quantities 
prescribed for the x direction only 

The Conservation of Momentum law, otherwise known as Newton’s Second Law, 
will now be applied to a fluid element. Newton’s Second Law states that the rate of 
change of momentum plus the net inflow of momentum within a control volume is 
equal to the net force on the control volume. For this derivation, we will only consider 
the x-direction. A similar analysis can easily be applied to the y and z directions. We 
will also only consider laminar flow. 
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Therefore, the momentum equation for the x-direction can be written as: 

d(M), 


W) out + (M) in + 


m dt 


-l ' 'I 


F] 


(B.l) 


Further expansion of the right hand side leads to the following: 

<W cv n 




m dt 


Jx 


surface pressure + body\x ^ ^ 


For our analysis, we will neglect any body forces which include terms such as forces due 
to gravity. It then follows after substituting the values shown in Figure B. 1 , the momen- 
tum equation applied to a fluid particle is: 

[(p vdxdz)u + 3(P v ^ xdzu dy + (pwdxdy)u + ^^ yu dz + 


, . , s d(pu)dydzu . 

(pudydz)u + — .. Ax 

ox 


- [(pudydz)u + (p vdxdz)u + (p wdxdy)u] + 


d(p«) _ 
dt 


d(o x )dydz d(x)dxdz 

o dydz + 3 dx + x xv dxdz H 37 dy + 


d(x xz )dxdy 

x dxdy + dz 


dz 


W.'V I V I -V 

dx x > dy 

- [o dydz + x dxdz + x dxdy] 


Further reducing we have: 

d(pu)dydzu ^ d(p v)dxdzu, d(p w)dxdyu , d(pu) 
d^c y+ dz dt 

d(o x )dydz d(x )dxdz d(x xz )dxdy 

3 dx + — -4r dy + 3 dz 

dx dy dz 

We can further reduce the previous equation and rearrange it to get: 

d(pw) d(p»~) d(pMW) d(pMv) _ *^xy ^xz 

dt + dx + dz + dy dx dy dz 

With the following relationship we can rearrange equation (B.5): 


(B.3) 


(B4) 


(B.5) 


°x = -P + 'x 


(B.6) 
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With this last relationship, we can now show that equation B.5 is the same as equation 

( 2 . 2 ): 


3(P«) , d(pu 2 + p) . 0(p«v) , 3(pMw) _ dX xx , . dX x z 71 

sr+ Tx 57 "57 "5z" 1 ' 

Furthermore, it can be similarly shown that for the y and z directions, the following two 
equations exist respectively: 



3(pw) 8(pMVy) 3(pVW) 3(pvr 2 + p) _ ^ T yz <Kz 

dt + dx + 3y dz dx dy dz 
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APPENDIX C 


Derivation of the Energy Equation 



Figure C. 1 Differential fluid elements with energy quantities 
prescribed for the x direction only 


The Energy Equation is derived from the first law of thermodynamics. Applying this 
law to a fluid element as shown above in Figure C.l„ the law states that the rate of 
change of energy inside the differential element is equal the net flux of heat into the ele- 
ment plus the rate of work done on element due to pressure and stress forces on the ele- 
ment surface. 
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In equation form this is simply: 


DE t 

~Dt 


= %e,~ W 


(C.l) 


Where E t is the total energy, q is the rate of volumetric heat addition into the element 


and finally, W is the rate of work done on the element. As with prior derivations, body 
forces such as those resulting from gravity will be neglected. Beginning with the left- 
hand-side of equation (C.l), the rate of change of energy within the fluid element is the 
substantial derivative of the total Energy, E t . Expanding this term out we get the follow- 


ing: 


DE dE t d(E,u ) d(E.v) d(E t w) 

+ - + J + 


(C.2) 


Dt dt 3 jc dy dz 
We will now move onto the first term on the right-hand-side of equation (C.l). This term 
represents the net flux of heat into the element. Referring to Figure C.l, we will add the 
appropriate terms to this element. For our discussion, the figure represent only the x- 
direction. A similar analysis can be done for the y and z directions. The net flux is then 
written as: 


( rXqjdydz , 

Mne, = q x dydz-[q x dydz + ^ dx J (C.3) 


Further reducing: 


dq x 

^ net ~ ~dx 


(C.4) 
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Adding in the y and z directions we finally get: 


M x M v M z 

[q] net = _ aF 


(C.5) 


The last term on the right-hand-side of equation (C.l) is the rate of work done on the dif- 
ferential element. Again, referring to Figure C.l, we will again make the appropriate 
substitutions. This term can be written as: 


W = 


updydz-{updydz ^ {UP l d x yd ± d*)] + 

dx j - X x .^udydz\ + 


x udydz + 


X udxdz + 

xy 


V 

d{x xx u)dydz 

dx 

d(x u)dxdz 


xudxdy + 


dy 

d(XM)dxdy 

dz 


dy^-x xy udxdz J + 
dz j - X zx udxdy j 


(C.6) 


Making further reductions in the above equation we get: 


W = - 


d(up) 


M* xx “) d(* X v“) M* zx u) 


xy 


zx 


(C.7) 


dx dx dy dz 

We will now add the y and z direction components to the above equation. The final form 
of the rate of work term is: 


3 (.up) , . ,){x *? u) 

w = --dr + sr- it & 

3 (VP) . , 9( V V) 

dy dx By dz 

3 (wp) S (V> . <X T yz w ) , d(x zz w > 

dz + 3x dy dz 
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We will now finally substitute equations (C.2), (C.5), (C.8) into equation (C.l) and get: 


dE t d(E t u ) d(E ( v ) d(E ( w) _ dq x dq y dq z 
dt + dx + dy + By dx By dz 

d(u P ) t t g(yo , a(T ^ M) 

Bx dx dy dz Q v 

_ B(vp) d(T xy v) d( T y V v ) d(T yz v) _ 

By dx dy dz 

d(w P ) . fly*') > , 3( \ z w) 

dz dx dy dz 

Rearranging this equation, we can now be shown that equations (2.11) and (C.10) are the 
same: 


dE t d((E t + p)u ) d((E f + p)v) d((E t + p)w) 
"5T + dx + dy + dz 


Bx 


+ 


3 ( M V + VT yy + M, V-gy) 
dy 


+ 


(C.10) 


+ VT vz + ^ ~ gP 

dz 


Page 66 



APPENDIX D 


Nondimensional Form of the Navier-Stokes Equations 


We will begin our analysis by first rewriting the governing equations. The following 
are the continuity, momentum, and energy equations respectively from appendices A, B, 
and C: 


dp djp u) tj(py) d(pw) = 
dt dx dy dz 


d(pu) , d{pu 2 + p) _ 3(pMv) , d(puw) _ ^ X xx , ^ X xy , ^ X x z ^ ^ 

dt dx dy dz dx dy dz 

d(pv) , djpuv) d(pv 2 + p) 3(pvw) ^ yz fn - 

at ox ay dz dx dy dz 

djpw) djpuw ) 3(pvw) 3(pw 2 + p) j^yz , n 

dt dx dy dz dx dy dz K 1 


dE d((E + p)u) d((E + p)v) d((E + p)w) 

— : ^ -j- -I = 

dt dx dy dz 

^ UX xx + VX xy + WX xz ~ q x'> d ( UX xy + VX yy + % ~ q y ) 

dx dy 

dz 


(D.5) 


We will now nondimensionalize these equations one at a time. 
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The nondimensional variables will be the same as those in equation (2.21). They are: 


u* = 


v* = 


w* = 


1 i i 


(D.6) 


The asterisk terms are nondimensional. They can be rewritten to solve for the dimen- 
sional variables. 


x = x*L 

y = y*L 
z = z*L 
t = t*(L/VJ 


U = M*V 00 

V = v*V oa 

W = H’* V 0 

ft = 


P = P*Poo 

P = P* V cope 

t = r*7\ 


(D.7) 


Substituting the appropriate values from equation (D.6) into equation (D.l) we have: 
3(p*pj f)(p*p M «*V M ) 3(p*P M V*VJ a(p*p oo w*v oo ) _ 

-f- “f“ V (^•^) 

d(t*(L/VJ) d(x*L) d(y*L) d(z*L) 


The left -hand-side of equation (D.8) can be factored further. 

Poo^qo / ^P*) 3(p*a*) + 3(p*v*) + 9(p*w*) 
L V dt* dx* 5y* dz* 


With further reduction, the resulting equation is: 


8(p*) ^ 3(p*M*) | 3(p*v*) f 9(p*w*) = Q 
dt* dx* 3y* dz* 
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Next, we will nondimensionalize the x direction of the momentum equations. Again, we 
will substitute the appropriate values from equations (D.7). We will first begin with the 
left-hand-side of equation (D.2): 


d(P*Poo w * V cJ d(P*Poo u * 2v2 °° + P* v JpJ 

a(rWVJ) + d(x*L) 

a (Poo M * V oo V * V oo) BiP^V^VJ 


(D-11) 


d(y*L) 


d(z*L) 


= LHS 


With further factorization, we get: 


2 2 

Poo^oo Vd(p*M*) + 3(p*w* + P*) + d(w*v*) + d(w*w*) 




3jc* 


dy 51 


3z* 


= LHS (D.12) 


Next, we will do the right-hand-side of equation (D.2). It follows then: 


RHS = + 


d(x*L) d(y*L) 3(z*L) 


(D.13) 


With further factorization, we get: 


RHS = |IY^ + ^ + dt *« 


dy* dz* 

The nondimensional viscous stress tensors follow as: 


(D.14) 


*xx = 3 ^. 


d(«*V ) ) d(w*V) 


d(x*L) d(y*L ) 3(z*L) 


(D.15) 


With further factorization, we get: 


_ 

1 XX 


dv* dw* 
3z* 


H 2 


(D.16) 
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Also, 


^ X y = 


( d(u*VJ d(v*Vj 
d(y*L) d(x*L) 


With further factorization, we get: 


r* = 

“ xy 


Ll V 

r oo oo 


M; 


3m* 3v* 

1" 

3y* 3** 


Also, 


T* = \L*)l e 


p(«'*V 00 ) 8(«*VJ) 

k dix*T) + d(z*L) 


With further factorization, we get: 


xz 


= (^HI^ + 19 


Substituting equations (D.16), (D.18), and (D.20) into equation (D.14) we get: 


RHS = , + ^ + 


3.x* 3 y* 3 z* 


p V 

If we divide both sides by _^_r_ we get: 


3(p*M*) + 3(p*M* + p *) + 3(m*v*) + 3(m*w*) 


3f* 


3jc* 


3/ 


3z* 


= LHS 


RHS = 


f L VlV ^-V- Y ^ t t xz 


vPooV^ ; 


3x* 3y* 3z* 


(D.17) 

(D.18) 

(D.19) 

(D.20) 

(D.21) 

(D.22) 

(D.23) 
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With further reduction and defining the Reynolds number as: 


re l = 


p V L 

• oo oo 

Hoo 


(D.24) 


The x direction momentum equation follows as: 


3(p*«*) d(p*w* + p*) + d(u*v*) + 3(w*w*) _ 

5/* 3x* 3y* 3z* 

1 vgg I dx **y I 

/?£jrA dx* dy* dz* 


(D.25) 


With a similar analysis, it can be shown that the y and z directions respectively are: 


d(p*v*) 3(p*wv) 3(p*v* + p*) 3(v*w*) _ 

dt* 3;c* 3y* dz* 


1 V dT *xy dT *vv dT * 


yy 


yz 


dx* 9y H 


dz* 


(D.26) 


9(p*w*) 9 (m*h’*) + d(v*w*) + 3(p*w* + p*) _ 

dt* dx* dy* dz* 


**xz , **yz , **zz 

dz* 


RE^Jy dx* dy H 


(D.27) 
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Lastly, we will non dimensionalize the energy equation. Again, we will substitute the 
appropriate terms from the list of equations in (D.6). 


dE* ( a((£*,t p *p„v 2 „)«*vj 

WVJ) + d(x*L) 


+ 


d((E*, + p* pjl)v»VJ d((E* l + P *p„V 2 ~)w*VJ 
3(y*L) d(z*L ) 

+ v* V^T* + w* V^T* - q* x ) 

: + 

3(jc*L) 

d(U* Vjt*^ + V* VjZ*yy + W* V oo X* yz ~ 

3 (y*L) 


3(m s| 


V T* 

V OO 1 XZ 


+ v*V oo T* yz + W 


V ~ x *z z-1* 


<Kz*L) 


Where 


E t = P [ e + 


2 2 2 

u +v +w 


Or in nondimensional form: 


E* t = p*p c 


f J 2 *2 2 *2 2 ^ 

^ . ,2 u* V oo + v* V oo + w* Vo 
e*V oo + 


V 


) 


(D.28) 


(D.29) 


(D.30) 
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Reducing, we get: 


2 2 2 

r->5t. 17 2 ^*f * , M* + V* + VC* 

E*t = (Poo V ~)P [ e + o 


(D.31) 


Also, from previous derivations, the nondimensional viscous stress tensors follow as: 


^oo^oo\2 „/ 0 3 m* dv* 3vv* 

T*_ = [—7— h\L* 2— - — - 




L )3 V 3z* 


(D.32) 


x yy ~ 


V~ v ~\2 3v* _ 3 m* _ 3 w* 

3 y* 3jc* 3z* 


H 2 


(D.33) 


r* _ 

zz 


J„dw* 3m* 3v* 
3z* 3 jc* 3y* 


M 2 


(D.34) 


r* = 
1 .ry 


Poo^ooA .fdu* 3v* 


MS* 


dy* 3jc* 


(D.35) 


r* = 
L xz 


u. V 

I oo c 


/ v3x* 3z* 


(D.36) 




u V 

* CO c 


r* - ■ ” | ^ 


: M 


3z* 3y* 


(D.37) 


Lastly, the heat flux term in the x direction is given by: 

- h— - 37 

q x~~ k dx~ (y - l)Pr3x 


(D.38) 
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In nondimensional form we get: 


YJll*|L*oo Rd(T*TJ 
q x ~ ”(Y- 1 )Pr 9(jt*L) 

Rearranging we get: 

l^ooY^T’oo dT* 

" * ■ -ITrTfTf 

It then follows the heat flux terms for the y and z directions are: 


(D.39) 


(D.40) 


q* 


y 


JT* 

LPr(y- 1 d y * 


i^y RT ^ a t* 

a* = U* 

H z LPr(y-ly dz* 


(D.41) 


(D.42) 


Next, we will substitute equations (D.31-D38) and equations (D.40-D.42) into equation 
(D.28). 
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With some factoring, this results in: 


(p.y~w*, poo^ 3 ” vtcg*, + p*)«*> i 

v 1 J'^ r+ l l J ar * 

fP«, V ' 3 ~V(( £ *( + P*)V*) ( Po= V ' 3 “y<( £ '*/ + P* )**'*) 

-Z - J a? + l~ J a? 

* v*T%y + gg*L> ^ I 

K l} J d** 

+ V***yy + + 

, L 2 J 3 >* 

+ v *'^* yz + tt ’*' r *ZZ> 9< ?*Z 

, z . 2 J 3z * dz * 


(D.43) 
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If we divide both sides by anc * using the definition of the Reynolds number, equa- 

tion (D.43) becomes: 

dE* t d((E* t + p*)u*) d((E* t + p*)v*) d((E* t + p*)w*)_ 
dt* dx* dy* dz* 

( 1 ^( M * T *xx +v * T *^ + vv * T:t: .rz-^V | 

\RE l j dx* 

(D.44) 

/ i \d(u*X* xy + v*T* yy + w*T* yz -q y *) 

{ rej) + 

f 1 \ 3(u*t* xz + v*X* yz +w*X* zz -q*) 

\RE l J d z * 


Where the heat flux terms are defined as: 

» _ H* ' )T * 

<7 x 2 * 

RE l M m (y-1 )Pr dx 

* M* 37^ 

q y RE L Mj(y-l)Pr d y* 


(D.45) 


(D.46) 


<?* 7 = 


(D.47) 


RE l M„ (Y-1 )Pr 


From equations (D.9), (D.24-D.26), and (D.42), we can place these in equation in vector 
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form. The resulting expression is the same as equation (2.35). 


| dE* ^ dF* ^dG* _ 1 ^\ + ^* v + 3G_ v 

ar* - + + + 3z* ~ ^ L V 3 jc* dy* dz* 


Where £?*>£*> F*, G* are: 


Q* 


p* 


p*w* 


p*V* 


p*w* 

p* «* 


p*w*^ + /?* 


P*M*V* 


p*W*U’* 

p* v* 

£*= 

p*a*v* 

F* — 

p* v *“ + p* 

G* = 

p* y* yy* 

p * w* 


p *u*w* 


p* y* vv* 


P*W*~ + p * 

£*, 


(£*, + /?*)v* 


(E*, + p*)v* 


(E* t + p*)w* 


And E* v , F* v , G* v are: 


E*„ = 


0 


xy 


u* X* xx + v*t* xy + w*X* xz - q’ x 


F* = 


G* v = 


0 


1 xy 

1 yy 

-r* 

1 yz 

M* X* xy + v*x* v> , + w*X* yz - q* y 

0 

1 X z 

T* 

i >>z 


u*X* xz + v*T* yz + w*T* zz - q* z 


(D.48) 


(D.49) 


(D.50) 
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APPENDIX E 


Transformation of the Navier-Stokes Equations 


First we will begin by solving equation (2.64). Repeated for convenience: 


10 0 0 

*>, ^ 

% i\ x n y 

Ci C, Cv C* 


10 0 0 

^ 

^ ^ y c 

z x z^ z^ z ; 


In order to solve for each matrix element on the left-hand-side, it will be necessary to take the 
inverse of the right-hand-side matrix. To expedite this calculation, a symbolic math program, 
Maple V was used. Figure E. 1 is the output from this Maple V session. 
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[ Si *T 

dutOa tty&dta - dsuCZMo st&LO* ittiOo dfdZ* to - i UP Z tta ttydUa | 

C1 — -. | 

[ -dtduat dycfTi tvSTnn . .ii&ou sfydTeta di/ST. - dydsau du/J 1 1 ) yl7.no - rtyd\ <m dtdZeia dtftXx - dytlnu dvtTs ttydZne * dfrdtau drd7.no dyrfT) -dtdTs didlna ♦ rfynTrta rfwQTi 

t ' * *r ” ■ ^ 

-JvtX\ iidZna * isttZna iuiXi OmUS dntZna - ixdZno drdX) ] 

_____ _ j 

f dxdioii dptX\ - itubui dydtto diAX> - dr&xu iviX\ ii&Zta - d*Bw dvltlo tofX . AzAau d.*JX. drd£i* drdUcu ixdHa ty£X\ -dydXi lUsiHo ♦ ifrfLOc dJdXi 

| %I * 

dxdXiiUd£ta rdtAEsadjjtXi dxiXs. dpUM . dtOOg d*Oi i j 


% I . dxdXi dydtio cUjtZna - dtdXi don.no sLjLtu - dydXi tUdlia dytZna » dydXi di wErta (UJ.na . djJXi dodlta dydlrta - djjSXs drdleta dydPJa 


Figure E. 1 Maple V session, calculation of metrics. 
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With the following generalized transformation, the governing equations will be transformed. 


* 

II 

Tv 

tTt 

(E.2) 

y = yCS.'n. 0 

(E.3) 

Z = z(& ri, 0 

(E.4) 

t = X 

(E.5) 


Using the chain rule of partial differentiation on these expressions, the following deriva- 
tives result: 


d_ 3|_3_ ari_a_ ag_3_ 

dx dxdt, + 3x3r| + dxdC, 


_3_ 3ti_3_ d^d_ 

dy 3y3% + 3y3r| + 3y3^ 

3 3£_3_ 3r) 3 3^ 3 

3 Z ~ a z a^ + 3z3Ti + 3za? 

3 _ 3^ 3_ dn 3 35 3 3 

3 ( 3j35 + 3(3ii 3(35 3t 


These derivatives will then have to be applied to the governing equations. For this anal- 
ysis, we will apply the derivatives separately to each of the nondimensional elements of 


equation (2.35). The continuity equation written in non dimensional form is: 


3p 3(pa) 3(pv) 3(pw) = Q 
3 1 dx 3 y 3 z 


(E.10) 
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Applying the transformation derivatives we have: 


cgap an ap a^ap ap 
aras ardn a/ac ax 

a^a(p u) aria(pM) a^a(pw) 
dx a^ + dx an + a* ac + 
a^ a(pv) ari a(pv) a^ a(pv) 
ay a^ a y ari ay a; 
a^ a(pw) an a(pw) a^ a(pw) 0 
dz dt, dz dr\ dz 


Rearranging and reducing: 


+ n 3p + r ^P + ^P + 

1 'dri at 

t 9(P«) . „ jj(gjO , r d(P«) , 

S a^ +T| * an - r ac 

t ^(pv) . ~ d(pv) , r 9(pv) , 

S' a^ S art S a^ 

C d(P*0 a(pw) r 3(pw) n 
S a^ S an S a; 


ap 

at 


, ap . a(pt<) , 

^ 3^ S y 

ap d(pji) 

ti.v- + n \ + n 

r an x an y 


a(pv) ? a(pw) 

a^ s a^ 


a(pv) + a(pw) 


an 'y an ' z an 
r ap r a(pn) . r acpv) r a(pw) _ 0 

Sf-jr T S r V. nr ' r V nr 


5 '3£ 


ac 


ac 




(E. 11) 


(El 2) 


(E.13) 


at 


a[p(^+^ M + ^v + ^w)] 

31 


atpOi, + V* + TV'V*')] 
an 

atpcS + ^M + C^v + c,^)] 

a? 


+ 

= o 


(E.14) 
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As mentioned in chapter two, the quantities in parenthesis of equation (E.14) are the con- 
travariant velocities. Restated, they are: 


U 

V 

w 


%, + 

ti,+iy + v + v 

+ Z x y + C/ + & 


(E.15) 


Rewriting equation (E.14): 


dp . atp U) . 8(pV) 9(pW) _ 0 

to % an a? 


(E.16) 


Next, we will transform the x direction momentum equation. In non-dimensional form, 
the x component of the momentum equation is: 

2 ' ' •*'- ^ ' 1 (to XX . dX xy^xz' 


djpu) 3(p U~_+_p) 3(pKV) 3(P«H>) _ 1 f^xx | ' 

dt dx dy dz R e i \ dx 


dy dz 


Applying the transformation derivatives to equation (E.17): 

3§3(£m) 3t| 3(pm) 3C 3(p«) 3(pu) 
dt dt, dt 3t| dt 3 C dx 

3^3(p u~ + p) , 3r|3(pM~ + p) , 3£3(p« + p) . 

rr T T ' 


dx 31; 


dx 3r| 


dx dt, 


3^ 3(pMv) 3r| 3(pMv) 3^ 3(p«v) 

3y 3^ + 3y 3t| dy dt, 

d£, d(puw) dr] d(puw) 3£ 3(p uw) _ 
dz dt, + dz 3q dz dt, 


1 1 

(dt^XX 


d^xx 

ReJ 

Hx 35 

+ 3x 3r) 

+ dr i) C 

1 | 

(S&x, 


. Z&xy 

Re L ' 

U y 3$ 

3 _y 3q 

+ 3y35 

1 , 

(d 

3n^ T Jt Z 

. 

Re L 

(dz 31; 

+ dz art 

dz 3; , 


(E.17) 


(E.18) 
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Rearranging and reducing: 


Further 



U d(p u 2 + p) . „ 3(p u+p) . r 8(p m 2 + p) . 

^ ac ^ an ^ a; 

e a(pwv) . „ acpMv) , r a(pMv) A 

q >' ac ^ an ^ ac 

e 3 (p«w) , „ a(pMw) , r a(pMH-) _ 

^ ac an ^ a; " 



d(pu) 

at 


t a(p«) . t a(p u+p ) , E a(puv) , p a(pMw) , 

^ a^ ^ a? + S’ ^ +< =z 3 E 


3^ 


3^ 




dCpM> . _ <Kpu2±P) . - 3(p«v) 3(pnw) 

"'“ail 11 * r)n +T1 >’ an ^ an 


an 


an 


an 


r 3(p«) . r a(p u+p ) . r a(pwv) . r a(pMw) _ 
^ ac ^ a^ ^ ac ^ ac 

1 ^ 3^ g 1*2- ax - 

' 1 " Sv lit Sz 


Re z 

_1_ 

Re, 


Re, 


3^ 


* ac 

3t 3t 

’vsrv^ 


3Z 3C 
z an 


3 t rr at rv „ at r , 

r — + C — s* + r « 

'sj: T 'sv T > 


3C 


ac ^ ac 


(E. 19) 


(E.20) 
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Also 


d(P“) . 

3t 

3[pw(^, + u + % y v + % z n) + k x p] 

3[p»Ql, + r|^ + ly + n^ + ^P] | 

dn 

dipu^+^u + ^.v + ^ + ^p] _ 

ac 

1 a ^ x ,, + fyty + ^ x «> 

Re L dt; 

1 3(V„ + V^ + ^ X «) . 

3^ 

1 a (^jr X x* + ^y^xy + ^>z^xz) 

Re L 3C, 

Once again, the contravariant velocities can be reduced: 

3(p«) 3(p uU + Z, x p) 3(p«V + r\ x p) 3(p uW + ^p) 
IT* 35 + 3n + a? 

1 3(^T„ + 5,t„ + 5 t T ja ) 

i 3 (n^„+V« + V«> 

<« t 34 

3£ 


(E.21) 


(E.22) 


Also, the viscous stress tensors are given by: 

_ 2 frydu dv 3w^ 

X xx - 3^1 Tx~Ty~Hz) 


2 (~dv 3 u dw\ 

x yy = 


(E.23) 

(E.24) 
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-2 (%— dw _ 3v 
T zz 3^V dz dx 9y 


(E.25) 


, du 3v 

x *y ~ ^ + Tx 


(E.26) 


xz 


fdw du 
~ Hfc + dz 


(E.27) 


fdv dw 
V ” + dy 


(E.28) 


These terms will also have to be transformed. Applying the transformation derivatives 
to equation (E.23), we get: 


X 


XX 


= U2m + 2^ + 


3T dxdt, dxdr | 


2 d^du] 

dxdO 


2 fd^dv dr) dv 

3^al + ayaTi dy^ 



(E.29) 


Rearranging: 



(E.30) 
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Transforming equation (E.24): 


x = u4% +1 ^ + m)- 

yy 3 P v i dyd^J 

2 fd^du dr\du 3C,aw3 

3^15x51 + afari + dxdO 

2 fdt,dw 3t| 

3 ,i v5ia? + 3z5n + 3iacJ 

(E.31) 

Rearranging: 


2 (~ e. dv t du c dw\ 

x yy ~ 3^1 2 ^yd^ x d^ z dO 

2 ( ~ dv du BvtA 

3^l 2 ^>’dn -T|jc 3Ti _Tlz 9riJ + 

2 f ^ dv - du * 3wA 

5^1 ^ydtr^dC^^J 

(E.32) 

Transforming equation (E.25): 


2 0 3t|3w _3^3 wA 

T =-3< 2 5^ + 2 fe5n + 2 5^J _ 

2 fd^du fr\du d^du]_ 

3^\d-*d^ + Bxdri + a*a^J 

2 fdt,dv 3r|9v 

3^a| + ayan + 9T3d 

(E.33) 

Rearranging: 


2 f. t 3w t du & 

T zz “ 3^l 2 ^a^ ^ ^aj 

2 (~ dw du BvA 

2 („ y 5vv r du r dv\ 

(E.34) 
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Transforming equation (E.26): 


fd^du dr\du d^du\ fd^dv dr| dv dl,dv\ 

“ ^5555 + 5 ?drj + 5 ja?J + H lac 5 | + a^ 5 n + 


Rearranging: 




du - du t 3v 3v 



Transforming equation (E.27): 


T 


xz 


= fl 


d^dw 3r| dw 

aIa? + a^ari + ^95J 


+ ^i 


d^dw drjdw d^diA 

aid| + fe arj + diW 


Rearranging: 


1 


xz 


(y dw dw y dw t du du y du\ 

+ n *3n + ^ + ^ + ^ + ^rd 


Transforming equation (E.2): 



d^dv dr\dv d^dv) f 

&al + 3laii + ai3CJ + ^ 


d^dw dt| dw d^dw A 

Tyli% + + TylX,) 


Rearranging: 


f d dv dv y dv t dw dw y dw 

V " ^d| + ^zdr j + ^dc, + a? + ^dri + ^ d£ 


A similar analysis can be done for the y and z component momentum equations. 


(E.35) 


(E.36) 


(E.37) 


(E.38) 


(E.39) 


(E.40) 
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The transformations for these equations are respectively: 


9fov) 3(P vU + Z, y p) 3(pvV + Ti yP ) 3(p vW + $ yP ) _ 

"aT + — 3^ + ^ + ac 

i + _ 

Re L 3^ 

1 ^,t cy + q v T v> , + ri z x vz ) 

Re L ' ^ + 

1 ^^A T jcy + ^>’ T >’>’ + 

Re L dt 

atowi 0(p wU + ^,p) 3(pwV + Tl,p) _ dip wWj^p) _ 

aT + aT^" + 5n * 

1 XZx'n + tyyz + tezz) 

Re L 

1 dtoS^ + Wy' + T)^) 

Re L 3^ 

1 d(^« + C y T y z + ^z) 

Re L dC 

The last equation is the energy equation. Restating the energy equation we have: 

dE f d((E f + p)u ) d{{E ( + p)v) d( (E t + p)w) _ 

~di + 5jc + 3y dz 

1 + vT xy + ~ t 1 3( “ X *y + vx >v + wX yz ~ SA + 

Re L dx Re L dy 

1 a (^z + vT vz + H;T zz-^> 

Re h Bz 


(E.41) 


(E.42) 


(E.43) 
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Beginning first with the left-hand-side of equation (E.43), the transformation begins as: 


dE t d E. dE t dE. 

^>tW + ^17 ZfW + 17 


% 


dt, '*dr\ ^a; ax 

d((E t + p)u) d((E t + p)u) ^ d((E ( + p)u) 


+ ^1 




x ' x arj ' a; 


a((£ f + p)v) a((E. + p)v) a((E £ + /?)v) 

^ 51 + ^ an + ^ 5? + 

a((£ r + /?)w) a((^ + /?)w) ^d((E t + p)w) 




arj 


dC 


= LHS 


(E.44) 


Rearranging: 

dE. 

a 7 + 

d{^ t E t + ((£, + p)^ x u) + ((£, + p) 6 y v) + {{E t + ( 

w~ + 

a(Ti t E + «£, + p) yO + ((£, + p)y) + ((£, + p)^)) 

arj 

3 ( 1 , £, + ((£,+ pK x u) + ((£,+ J>) 1 /) + ((£, + PK Z W)) _ ( u< , 

51 

Further rearranging: 

dE { 

a 7 + 

3 ( 1 , E, + (E, + ?)(!,« + 5 / + l z w)) 

31 

3 (T|,E, + (E, + P)(n X « + * 1 / + T 1 ; W )) 

arj 

W,E, + (E t + p)& x u + t; y v + ty)) ruc 

31 


(E.45) 


(E.46) 
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Adding and Subtracting like terms to help in reduction: 


d(^ t E t + {% t p - % t p) + (E f + p){% x u + l y v + $ z w)) ( 
d(r\ t E t + (t \ t p -i\ t p) + (E ( + p)(r\ x u + ^1/ + ^)) 


(E.47) 


d(^ t E t + (£,/>- ^ t p) + (E { + p)(^ x u + £ y v + ^w)) _ 


= LHS 


Rearranging terms to help to combine terms defined as contravariant velocities we have: 


9((£, + p){% f + %u + % y v + ^ z w)-^ t p) 

w~ H 

d((E t + P){ T\ t + T| U + T|/ + \W)-T\ t p) 

an 

d((£, + />)(£, + + CyV + ^wK^) 

"aT” 


= LHS 


Final reduction yields: 

dE t dg^ + pVJ-^p) 

dr + 3£, 

d((E t + p)V-x] t p) _ d((E t + p)W-^p) 

3rj + 3£ 


= LHS 
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Continuing on with the right-hand-side, we have: 


* 35 +11j: an 

% uX xx + vx xy + WX xz “ . B a(uT xy ^ 

" + Sv 


ac 


a*; 


RHS(Re L )= 

^ d( MT x x + vT jp + wT sz~ gJ . „ ^“ T xx + vX xy_ +WT «“^ 
C ; 

T1 

5 

c 


a ^y + VT^ + WX^-^) ^V+VTyy + WX^-gy) ^ 


3ti 


a; 


d( U X xz + vly Z + WT, Z - q z ) d(ux xz + VT VZ + wx zz - q z ) 


+ % 


xz yz 


z d % ' ’ z 

^ uX xz + vX yz^\z~l±l 


an 


ac 


(E.50) 


Combining terms: 


RHS(Re L )= 

d &x<- ux xx + vX *y + wt xz-‘<x>'> , . 

si 35 

a ^ z (“ x x Z + VT v z + wx zz-1z>'> 

35 

3(T) v («T tx + VT JV + »t n - <?*)) 1 3(T| v (Mt xy + VT y) , + w^yz ~ fpj + (E.51) 

dll 3n 

3 ( t ' z ( mt x Z + vx yz + wx zz ~ ^ Z >> . 

9t| 

3(5,(«t„ + VT xy + »^-^)) _ 3(5 v (^ y + ^vv + ^v Z -‘?y)) , 

35 35 

a (^< HX x Z ' | - VT yz + K ’ T ZZ - g z )) 

35 

At this point we have already transformed the viscous shear stress terms. The only other 
terms that contain derivatives are the heat flux terms. 
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Applying the transformation to the x-direction heat flux term, we get. 


M- 


Q x ~ 2 

(y - 1 )M oo Re L Pr 


(d^dT dr\dT d^dT 
+ dx + dx d^ 


Reducing to its final form: 




(y- 1 )M oo Re L Pr 


' dT dT r dT 


Similarly, for the y and z directions respectively are: 

M- 


q y = 


= 


(Y- 1 )M ooRe L Pr 




t dT dT r dT 


e dT dT y dT 


(E.52) 


(E.53) 


(E.54) 


(E.55) 


(Y- 1)M oo Re L Pr 

Finally we can put equations (E.16), (E.22), (E.41), (E.42), (E.49), and (E.51) back into vec- 



Where 


Q = i 


F = 


-1 

pu 


puU + ^ x p 

E = -J 

pvU + t, y p 

? 

P wU + Z, z p 

- 

( E f + p)U - ^ t p 


pV 


pW 

puV + r] x p 

o- l l 

puW + ^ x p 

pvV + r\ y p 

pvW + ^ y p 

pwV + TJ z p 


pwW + C, z p 

(E t + p)V-r\ t p\ 


(E t + p)W-^ t p 


The viscous flux terms E v , F v , G v are: 


E v = J«A + W 

= j(.r\ x E v + n y f v + \G V ) 

o v = j«A + ?A + w 

Here, both sides of the equation (E.54) have been divided by the Jacobian. 


(E.57) 


(E.58) 
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APPENDIX F 


Derivation of the Thin Layer Approximation 


Let us first begin with equation (E.22), the x-direction momentum equation. Written again: 




dr 

ft 

- + 

1 

+ lAy + 

^z X xz) m 

Re L 

a^ 


i 

+ V*y + 

■V«> 

Re L 

a$ 


1 

+ ^v + 

C z t«) 

Re L 

a^ 



*1 




(F.l) 


Lets only consider the right-hand-side of this equation, the viscous flux terms: 


i 3(^t rr + § v t J(V + ^x xz ) t + + 

RHS = rT l — dT + 5? 

1 ^x X xx + ^>' T *>’ + ^>Z X xz) 

Re, 3; 


(F.2) 


As stated in chapter two, the thin layer approximation reduces down to terms that contain deriva- 
tives only taken with respect to Zeta. Reducing equation (F.2) we have: 

1 + ^ y x *y + ^z x xz) 


RHS = 


Re, 


ft 


(F.3) 


Expanding this out we have: 




(^ + Ti z w n + C z ^))J + 

+ + + + + + 

+ ly*, + C. ," t + + v, + 


(F.4) 
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Further reduction yields: 


rhs = 


(F.5) 


+ C^)) 


/ac 


Expanding out: 


RHS = 


+ C,yU^ + CjtCyV^H + + CxCjK^M- J 


Re, 


(F.6) 


3C 


Rearranging: 


R//S = 


n( 5 + 4, + 4,1" Wj"';!! + CxC, v c^ + 4,4, : B !^ 


(F.7) 


Re, 




Further rearranging: 


R//S = 


H(44 + 4, + 4;)“; + ^‘4,4, ‘V + ^4,4;“'^ + 5^4>; 


(F.8) 


Re, 


ac 


Finally we have: 


RtfS = 


\l(£ + £ + £)«<; + \^x(^ + 5/C + ^ w 0] 


(F.9) 


Re, 


ac 


A similar analysis can be done for the y and z components of the momentum equations. 
There results respectively are: 


RHS = 


v&l + C>. + & x x + + + 

K 


Re, 


(F. 1 0) 
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RHS = 


\X(£ + Cy + + + V? + 


(F- 11) 


Re, 




Lastly, we will reduce the energy equation. Rewriting equation (E.49): 


RHS(Re L )= 

"V + VX xv + WX xz -^K 8 ^v( UT v> + VX vv + H ' T K~ g >' )) + 

51 31 

a (M“ x « + vT yz + wl zz-‘>z )) . 

31 

3(ll t (Ht xr + VT lyi + - <? ,)) 3(T| y (»tj 0 , + V ' I yy + tl ' X yz + (F. 12) 

3rj 5ii 

3(T1 z («^ z+ VX y .- | -H’t zz -< ?z ) ) | 

3ti 

ad^x^ + vx^ + wx^-y) _ a( y«x^ + vx >iy + »'x yz --7 y )) i 

— 5i 31 

a(l z («x :cz + v Xy Z + wx zz -<?.)) 

31 

Reducing yields: 


RHS(Re L )= 

a<ly(" x vy + vx IV + H ' x «-9y) ) , a <l y ( “ x yy + t ’ x v.v + H ’V-gy>L (F ,13) 

—— at 31 

3 ( 1 z (mX xz + vx yz + w x zz - ^ z )) 

31 
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Let us first reduce each of the shear stress terms and the heat flux terms. 



2 f _ t 3m t 3v e dw\ 

X xx - ^ad 

2 ( _ 3m 3v 3wA 

3^ 2 ^^ _ ^a^ _ ^3nJ 

2 du «, 3v ~ dw\ 

3^l^a? _ ^ae, ^a;J 

(FI 4) 

Reducing: 

2 3m r 3v r dw\ 

X xx - 3**v ^a; Sac, ^acJ 

(F.15) 

Similarly: 

2 3v r 3 m r 3wA 

(F.16) 


2 f^y 3w y du y 

*= = 

(FI 7) 


/V 3m j- 3v"\ 

X xy - K^d? + ^3?J 

(F. 1 8) 


+ C lS r J 

(F.19) 


3v r 3wA 

x yz ~ ^ac + ^aj 

(F.20) 

The heat flux terms can be reduced to the following: 



n ( r dT\ 

(7-l)M 2 oo/?e L Pr 

(F.21) 
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Similarly for the y and z directions: 


(r dH 

qy (y- \)M 2 ooRe L PA 

= ]i (y ?!) 

( y - 1 )M 2 ooRe L Pr K Z 

Substituting the shear stress and heat flux terms into equation (F.13), we have: 


RHS(Re L )~ 

( 

c 


^ ac) + + ^ac) + "'K^ac + ^ac) + 


> - . MD 




(y- l)Af ~Re L Pr 


)) 


(y-\)M 2 ~>Re L Pr^ } ^ 


M))) 


ac 



( r dw _ au3 av _ aw 3 . 2 ,j. r aw r a« 3v3 
"H^ac + ^acJ + V> v*dc + ^a J + l Vr^ac ^ac ^ v aJ ‘ 


(Y - l)M\>Re L Pr 


(<«]) 


ac 


Expanding: 


RHS{Re L )= 

A r 2 a« 2 , a»’ 2 r r 3w' „ r r 3 h r <*“ iL 

^"^ac _ 3 K "^^ v ac _ 3^“^^ v ac v *^ v dc ^ac ^ac ^^ac 


>ac 


>ac *ac 


(Y-l )M'„Re L Pr Ul = / 


r 2ar 

^ac 


ac 


~20m *. w ov h yiov * r r ou r r uvv . „r r WK , ... r^ urv 

unCvSr + “^Aar + + "i^Agr + ^af + 


3v 4 k 20v 2 


3m 2 


3vf 


3v 


K 23w 




ac 


Z 3C 


^C 


H r 23r 

' ^ (y- \)M 2 ~Re L Pr } ^J 


y> u 3>V rZOK ylOV r Y ow ** t*ZC7W f UK ^ 


yldu j,23v 


3w 4 r 23w 2 


3m 2 


dv 


JJL ^23T 

(Y-l )M Z ~Re L Pr '^J 


ac 


(F.22) 


(F.23) 


(F.24) 


(F.25) 
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Rearranging: 


RHS(Re L )= 

-Jy 2f fd M 2 3v“ 3 h'“\ 

V'Hac + 8t 3t J 

+ 0t 3C. 

1 + 1 

(K'^) _ 3 R “^''3C + V ^'3C + W ^ x ^>i) 

+ 

(F.26) 

+ 

) 

$[nl 

fdu 2 8v 2 8 w 2 > 

IdC 3^ + 3^J 

dri 

l + a 8?J 

+ 



I'Sm’’ 3v 2 3h' 2 ' 

l8C' t "8C 3C ) 

1 37' 

l + “^_ 

+ 

3? 

fly 2 3h' 2 '\ 2 r , 3« 2 y y 3v , y 3w - y 3w 

U^ir rr ^-^-3 + u ^k + 


where a is: 


(y- \)M 2 -Re L Pr 

Reducing even further we finally, with further reduction, we have: 


(F.27) 


RHS(Re L )= 

We can now put equations (F.9) , (F. 1 0) , (F. 1 1) , and (F.28) into vector form. One will notice 
there are no viscous contributions from the continuity equation, and therefore the first 
element of the viscous flux vector will be zero. We can now show the thin layer viscous 


(£ + Cy + + % + If) + «|] + + C 'l + + C ' 1 ’ + ^ | (F.28) 


x, 


flux vector S as: 
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0 

+ Cy + + C v v ; + 

V-(C K + Cy + Q V $ + + Cy V S + 

H( C y + ^y + | ( £,«£ + 
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